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FORWARD 


Aiming  at  understanding  and  simulating  mixed-type  flows  of 
interest  in  applications  to  transonic  aerodynamics  and  plasmadynam- 
ics,  mathematical  and  numerical  analyses  for  relevant  initial 
boundary  value  problems  and  finite  element  discretizations  are 
described  in  the  following. 

The  report  is  divided  into  two  parts: 

I.  Various  finite  elements  formulations  and  discretizations  for 
initial/boundary  value  problems  for  the  Tricomi  equation  in  the 
hyperbolic  domain,  work  carried  out  by  Mrs.  Frieda  Loinger.  This 
is  a continuation  to  Sara  Yaniv's  work  reported  last  year,  and 
the  schemes  developed  are  tried  on  the  same  sample  problems. 


II.  An  attempt  to  extend  "hyperbolic"  methods  of  analysis  to 
mixed  type  flows  with  hyperbolic-elliptic  shocks  and  ensuing  en- 
tropy conditions  for  stable,  accurate  numerical  schemes,  done  by 
Prof.  Michael  Mock  with  an  appendix  by  Lustman  and  Geffen  on  the 
impossibility  of  elliptic-elliptic  shocks  (a  fact  well-known  to 
any  aerodynamicist  and  rather  obvious  on  physical  grounds.) 

The  question  of  entropy  functions  and  inequalities  is  inher- 
ently related  to  the  well  posedness  of  the  nondissipative  formu- 
lation and  in  turn,  we  feel,  to  the  character  and  behaviour  of  the 
different  variational  functionals.  A further  clarification  along 
this  way  may  be  essential  for  the  ability  to  correctly  select 
the  physically  relevant  fields  from  the  multitude  of  mathematic- 
ally possible  ones  (e.g.  with  expansion  shocks)  by  full-proof 
a-priori  conditions  in  the  simulation  scheme. 
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INTRODUCTION 

A finite-element  formulation  of  Tricomi's  problem  using 
bilinear  isoparametric  characteristic  rectangular  elements  in  the 
hyperbolic  domain  is  given  in  [l]  with  two  sample  problems.  In 
this  report  discrit izations  for  linear  triangular  elements  are 
given,  with  calculations  for  the  hyperbolic  part  of  the  sample  prob- 
lems mentioned  above.  The  following  discretizations  are  given  in 
detail : 

I.  Global  application  of  the  variational  conditions  for 

a)  Cauchy  b)  Goursat  conditions  r^s,ulting  in  an  implicit  scheme 
for  the  whole  hyperbolic  triangle. 

II.  Local  application  of  the  variational  conditions  on  each 
element  for  a)  Cauchy  b)  Goursat  boundaries  resulting  in  an 
implicit  scheme  (or  marching  along  characteristics). 

The  order  of  the  approximations  is: 

2 

0(h)  at  internal  points  for  all  implicit  schemes,  0(h)  near 

2/3  • • 

external  boundaries  and  0(h  ) near  the  parabolic  line  y=0. 

For  the  explicit  (marching)  schemes,  only  0(h)  is  achieved.  The 
isoparametric  schemes  become  inconsistent  near  the  parabolic  line, 
due  to  the  singular  behavior  (coalescence)  of  the  characteristics 
there;  hence  the  isoparametric  triangles  are  changed  to  regular 
(x,y)  triangles  for  which  consistency  is  restored.  It  also 

facilitates  further  local  refinements,  (.possibly  using  the 

multi-grid  method ), where  they  may  be  needed. 

An  appropriate  combination  of  linear  triangles  near  the  para- 


— » 
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bolic  line  and  isoparametric  characteristic  triangles  away  from  it 
give  the  most  accurate  results,  with  the  marching  scheme  being  most 
efficient  and  easiest  to  apply.  Schemes  using  Goursat's  conditions 
have  been  tried.  Convergence  has  been  demonstrated  for  implicit 
schemes  [1],  but  not  for  explicit  ones.  Details  will  be  reported. 


Reference 

[l]  Yaniv,  S.,  Variational  Formulation  for  Formally  Symmetric  Prob- 

( 

• lems  and  its  Application  to  the  Tricomi  Equation,  Final  Scientific 

Report,  Grant  AFOSR  73-256,  July  1977,  Part  I,  also  Ph.D.  Thesis, 
Tel-Aviv  University,  Israel,  1978. 
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!•  The  hyperbolic  problem 

Given  Tricomi’s  equation: 


yep  ~ m =0 
J XX  ^yy 


(1) 


or,  equivalently,  the  first  order  syste 


em 


where : 


yUX  - Vy  = 0 


u - v =0 
y x 


( u , V ) = (tp  ,tp  ), 

x y 


(la) 


(lb) 


we  look  for  a solution  in  the  hyperbolic  region  D,  bounded 
by  the  parabolic  line  y = 0 and  the  two  characteristics 
and  f2  where  (figure  1): 


and 


ri : -1 


r2:  +1 


= n = 


= £ = 


X - !v3/2 

x 3y 


x + 3 y 


2. .3/2 


■lCx<0 


0<x*l , 


We  assume  that  the  problem  is  well-posed,  with  either 
Cauchy  (i)  or  Goursat  (ii)  conditions  given: 

(i)  ip(x,0)  = f 1 ( x ) <py(x,0)  = f 2Cx) 
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or  equivalently  u(x,0)  = f^(x)  v(x,0)  = f 2 ^ x ) 

given  on  the  parabolic  line. 

(ii)  ip(x,y(x))  = f(x)  , or  equivalently  u/y+v  = ^Cy)  given  on  the 
characteristic  I,  Cy  = [j(x+1)]^/j) 

! 

and:  <p(x,0)  = i^(x)  or  equivalently  u(x,Q)  = f^Cx)  on  the 

parabolic  line. 


a B 

(u,v)(x,o) 


Figure  I a '•  Cauchy  data. 


ix/y+v^y)  — 


Figure  lb  : Oouruat  prob. 


1 


2 . 

(2) 


(2a) 

(2b) 


j»  (2c) 


T 


Variational  formulation  in  the  hyperbolic  domain 
Define  the  functional 


J (u , v , A ) 


[yu^-v^+A (u  -v  ) ] 
J y x 

D 


dxdy 


+ J A(udx+vdy) 

r 

m 

where  D is  the  hyperbolic  domain  and  f includes  all 

r m 

surfaces  where  no  boundary  conditions  are  given:  T = r-,ur„ 

m 1 2 

for  the  Cauchy  problem  and  rm  = ^ for  the  Goursat  problem; 

A is  kept  fixed  on  F . 

y m 

J(u,v,A)  becomes  stationary  when  (u,v,A)  satisfy  the 
equations : 


2yu  = A 


2v  = A 


yu  = v 
x y 


u - v =0 

y x 


which  are  equivalent  to  equations  (1)  for  A(x,y)  C C‘ 


— * 
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3 . Finite  Element  Schemes 

The  hyperbolic  region  is  divided  into  triangular  isoparametric 
elements,  (figure  1).  In  every  element  we  take  trial  functions 
which  are  linear  in  £ and  n. 
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A . Interior  points 

Minimizing  the  functional  for  interior  points  we  get  the  foll- 
owing schemes: 


3 = 1 , . . . , 2N-2 

i = -N+l , . . . ,N-1- j 


3u  . . 
i >1 


J (u , v , A ) = 0 -* 


,3,  .4/3  12r,  . w1.7/3,l<  ...  .7/3.3  .10/3  3 ..Al.10/3. 

(_h)  -7t(ui>j+1+u..lj.tl)(7]  +2(=l  + 1)  +T^>  ~T0^  + 1)  ) 


+ (u . . . +u  . , . . 

1-1,]  i+l,] 


w .7/3  9.13/3  9 ,.  . .13/3  J ,....13/3. 

)(~J  -653  +TT0(;|‘1)  +li0(;i+1)  ) 


130 


+ u.i.(-237/3+|(3+l)10/3-f(3-l)10/3 


+ 18,13/3  9(.  )13/3__9(.  }13/3) 

65:  6 5 1 65^1  ’ 


65 


(3a) 


*.  » • . 
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(Ai-l,  j+l~Ai,  j-l)  + (Xi,j  + l-Ai+i,j-l) 


Ihe  scheme  (3a)  approximates  the  equation  2yu  = A at  the 
mesh  point  ^),  and  the  accuracy  is: 


2yu-Ay  = 0(h) 


J (u  , v , A ) = 0 


* <vio-i*vi*i,j-i)(<i-1)5/3U5/34<i-u8/3458/3) 

(3b) 

= FttAifj+l“Ai-l,j  + i)(35/3  + (j  + l>5/3438/34(3+1>8/3> 

5/3 

+ (Xi+l,j"Ai-l,j)(_2i  ^(j  + D8/34(j-l)8/3) 

+ (Ai+i,j-i-Ai,j-i)((j-1>5/3+:5/34(3"1>8/3-|38/3):1 


11 


The  scheme  (3b)  approximates  the  equation  2v  = Ax  at  the 

mesh  point  (x.  . ,y.  . ) and  the  accuracy  is: 
i » 3 i » 3 


2v-Ax  = 0(h2) 


3*i,j 


J(u,v,A)  = 0 


,3.2/3 

12  S'* 


*-rT73“CCvi,j  + l“Vi-,i,j  + i)(35/3+(j+l)5/3+|j8/34c3+l)8/3) 


+ (v 


i+i,j-vi-i,j)(-2^5/3+l(3+i)8/34^-i)8/3) 


(3c) 


+ ( V • 


5 / 3 , , . . .5/3.3,.  , .8/3  3.8/3. 


i + 1 


,j-i_vi,j-i)(j  +Cj-1>  Yi-D'  -p  /a)] 


= (u.  1 -,-,-u.  . .)  + (u.  . , . - u . . , . -.  ) 

1-1> 1 + 1 1,3-1  1,3  + 1 1+1, 3-1 


The  scheme  (3c)  approximates  the  equation  u^  = v^  at  the 

mesh  point  (x.  . ,y.  .)  and  the  accuracy  is; 

1 5 J 1 > J 


Uy  - vx  = 0 (h2 ) 


Remark : 

For  an  interior  point  (x,y)  the  accuracy  of  the  equations 

(2a)-(2c)  is  0(h  ) only  if  we  assume  that  the  value  of  j 

is  sufficiently  large.  This  can  always  be  achieved  by  taking 

3 2/3 

smaller  values  of  h,  because  y = (^-hj ) so  j increases 

if  h decreases  for  a fixed  interior  point. 


-C|  u> 


I 
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B . Boundary  Points  for  the  Cauchy  problem 


In  this  case  I’m  = and  we  have  the  following 

functional 


2 2 

•J^  (u,v,A)  = [yu  -v  +ACu  -v^)]dxdy  + 


y x 


A (udx+vdy ) 


riUr2 


mesh  points  on  : i = -N 

j = 1,. . . , 2N-1 
n = -1 


h)4/3  1 8 r ,1.7/3.3.10/3.  3...  n10/3,  27.13/3  27  , ..,.13/3. 

h)  "T^-N.j  + l^3  +5^  +T0(:I+1)  +130J  "130(:1  + 1)  > 


+ u 


-N+l,j 


(.j7/3.5|j13/3ti|_(..i)13/3tT97r<.tl)13/3) 


130 


+ ,,  , -7/3  3,.  , x 10/  3 6.10/3.  9-13/3  18,.  ,.13/3.  9 13/ 3 . 

+ U-N , j (~^  ~P  +653  -65(J-1)  +65(^  + l>  > 


— » 


r^i  t— > 


t 


,17/31,.  .,7/3.  3 ,.  .,10/3  3.10/3, 

+ U-N  + 1 , j - i ? +nr(3-1)  -t*>  ) 


,1,.  . .7/3  3/  . .,10/3.  3.10/3,  27  ,.  .,  13/3  27.13/3.-, 

+ u-N,j-i(2<:]-1)  V]-1}  +To:  +TW(^-1)  -TltfJ  )] 


- A W fA  A-N,j-l*A-N+l,j-l. 

-N,j+1  -N+l.j-l'  -N,j  2 )- 


The  accuracy  is 


Tv  . Jl(u’v’A)  = 0 
__  "naJ 


2yu-A  = 0(h) 

y 


,.5/3. 3. 8/3. 3#.  ..,8/3.27. 11/3  27, ...  11/3, 
V-N,j+l(j  +2:  +4C^+1>  +44^  "44C^+1  } 


, - . 5/ 3 9.11/  3 9,.  , 11/3  9 ,.  ..,  11/3  , 

+ V-N  + l,j("23  "22^  +44(3-D  +44^  + 1>  > 


, ..5/  3 3 / . .,  8/  3 -.8/  3 9 ,.  .,  11/3.  9.11/  3.  9 ,...  . 11/  3 . 

V-N,j("2j  ~33  “TT^"15  22:  +27(3+1)  ) 


-N+l , j-1 


, , • .,5/3, .5/33,.  . ,8/3  3.8/3, 
((:-l)  +3  +jj-(]-l)  “4  3 ) 


...  1 ,5/3  3 , . . ,8/3  3.8/3  27 , . .,11/3  27-11/3. 
-N,j-l((:l"1)  +7(3-l>  +Tl  +inr(3-1)  -W3  } 


.-A_n  . )(-2j5/3+|(j+l)8/3-|(j-l)8/3) 


~N+1 , j -N , j 


, ,,.5/3,,.  ,,5/3.3,.  ,,8/3  3.8/3, 


V 
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l*  • Linear1  Triangular  Elements  near  the  Parabolic  Line 

To  avoid  the  inconsistency  near  the  parabolic  line 
for  small  values  of  j,  (near  j = 0),  we  must  change 
the  trial  functions  near  the  parabolic  line  for  the  same  mesh 
points 


i 


x 


( i+^-)h 


1 


- r3h^2/3 

- (4hl) 


and  take  ordinary  triangular  elements  with  linear  variation  in 
x and  y,  and  linear  approximations  to  and  r2  near  the 

parabolic  line. 

Let  the  approximated  domain  be  D,  and  the  approximated 
characteristics,  and  f2  . 

For  j < 1q  the  trial  functions  will  be  linear  in  x,y. 
For  j > jQ  the  trial  functions  will  be  linear  in  £,n. 

The  variational  formulation  is: 


J (u , v , X ) 


2 2 

[yu  -v  +A(u  -v  ) ]dxdy  + 

y x 


A(udx+vdt ) 


For  j ^ Jq  + 1 all  the  elements  belonging  to  the  mesh  point 
(i,j)  are  isoparametric,  so  the  schemes  are  the  same  as  in 
chapter  3. 
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For  ] = Jq  the  trial  functions  are  linear  in  £,n  for 

linear  in  x,y  for 


For  j < jjj  - 1 all  the  elements  belonging  to  (i,j)  are 
ordinary  triangles. 


* 

1 


A)  Interior  Points 


».(Ih)^3[(u1>.+1*ui.lj.tl).I|(j7^(jtl)7/3t|j10/3.|(jtl)10/3) 


+ (u.A.  . +u 


1,  1.4/3  1,.  ..4/3  1.2/3,.  ,.2/3  3.7/3 


i + 1 


. . \ -I-  / X • “T  / O X ✓ • t \ *T  / O X • / O / . - v £ / O O • 

,3  Ui-l,j)  2(-r]  “TC3  (3-!>  "73 


35 


9.10/3  27.13/3.  27/..-.  .13/3. 

J 455J  455  3 ' 


k\ 


I 
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+ ui,j(-^4/3-il^-1)4/3-^2/3-^7/3 

+ 4(j  + l)^/3+  27^3/3,  27(j+1)13/3) 


+ (u. 


1 ,,4/3  , x 4 / 3 . 


( Ai-1 , j + l“Ai  , j-l)  + ( Ai,j  + l_Ai  + l , j-l) 


(u  ,v  , A ) = 0 


(vi,j  + l + Vi-l,j  + lHlj5/3+f<j+1)5/3  + 2fj8/3(j+:L)8/3) 


+ (v.  n . + v .)(-Ij2/3_- i(i_i)2/3_2i5/3__9i8/3__^7ill/3 

1 + 1,1  i-1,1  4J  12  : . P 20:  220] 


27. ....11/3. 
+ 220(^+1)  5 


r , 5.2/3  lf.  , >2/3  6.5/3  9,.+1.8/3  27,  . . ..  .11/3  . 27.11/3. 

”(;]-1)  -5:  +I0(3  + 1)  -ITO(:)  + 1)  +U03  > 


+ (v.  .+v. 

i 


,i+vi+i,i-i)(^2/34^-^  /3) 


- 1rn  x w3.S/3  3f.  ,5/3  9.8/3  9 ,. .,.  8/3. 

" TT[Ui,j  + l-Ai-ilj  + i)(53  V3  + 1)  +253  '20<3+1)  > 


+ (A. 


-A  ) )2/3  3.5/3  _9  • 8/3  9..  .8/3. 

,j  Ai-l,jM  33  r3  1)  P "W3  +40(3  X)  ) 


— 9 


^i+l,3-l~xi,j-l 


)(lj2/3_l(j_1)2/3)] 


- ( u , v , A ) = 0 

d 


2/3 


' . 3,.  .5/3  ,.  . . ,5/33.8/3  3,  .,,,8/3. 


i + 1 


w 1-2/3  1,.  2 / 3 3.5/3  9-8/3.  9/..n%8/3»  . 

,j“vi-l,3)(TO3  T:  “TO3  +TO(:  + 1)  5 (lc) 


W1  , . 2/3  , . u2/3n 

i+l,j-l-vi,j-l)(r(^  -(3-1)  )] 


• i • -,)  + (u.  ..,-u..,  . ,) 

1-1,3  + ] i,3-l  i,3  + l i + l, 3-1 


4/3 

— [u.  • ((3+l)2/3-(3-l)2/3)(8Cj+l)2/3+8Cj-l)2/3+14j2/3) 

1 , 3 

. . )((j+l)2/3-(j-l)2/3)((3  + l)2/3  + C3-l)2/3  + 3j2/3)  <2a) 


, + 1 


I 
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+ (ui,juVilju)*5((3  + 1)4/3-54/3) 


+ (ui  + 1,j-]+ui,j-i)5(34/3-^-D4/3)] 


- (A-  — A • • *i  ) + ( A • * n * 1 1 ) 

l-l, j+1  i,]-l  1,3+1  l + l, 3-1 


( u , v , A ) = 0 


, . w .2/3  , . , .2/3. 

+ ( v . . , +v . . , . , ) ( 3 -(3-I)  ) 

1,3-1  l + l, 3-I  J J 


RC(xi,j+i-Ai-i,3+i)((^+1)2/3-j2/3) 


+ (Ai+l,3"Ai-l,3)((j+1)2/3"(j"1)2/3) 


+ (A  . 


2/3  ,,2/3. 


i + 1 , i 


,j-l“Ai.j-l,<j  "<3-1)  >] 


^ ( u , v , A ) = 0 ■* 

*-  >3 

,3. 2/3 

2 V r , w,..,,2/3  . 2/3. 

r hi/3  C(vi,ju*vi-i,j*i><<3tl>  > 


v 
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b 


i 

' 


* (vi+i,j-vi_i,j)((3+1>2/3-(:-i)2/3) 

+ (vi  + l,i-l-Vi,j-l)(32/3-(j-D2/3)] 

= Cui-l,j+l“ui,j-l)  + (ui,j+l"vi+l,j-l) 


(2c) 


Taylor  expansion  gives  that  the  order  of  the  equation  is: 

0(h2)  for  j * - 1 

0 (h)  for  j = jQ 

lor  sufficiently  large  values  of  j.  For  small  values  of  j the 
accuracy  of  the  equations  can  be  less. 
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J 


, 1.4/3  1,.  ,.4/3.  1,2/3/.  , . 2/ 3 . 


(A-N,j+l“A-N+l,j-l)+(A-N,j 


,18.5/3  27  8/3  27  ,.  , 8/3,243.11/3  243  , . . 11/3. 

J +—J  +To(^  + 1)  +Ho3  -TTo(J  + 1)  ) 


-Pico 
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+ u-N,j-1(2i4/3-3(3-l)4/3+j2/3(j-l)2/3)] 


- /a  . w/a  _A-N,i-l+A-N+l,i-ls 

- ^A-N,j  + l A-N+l,j-l)  + (A-N,j  2 ) 


— Ji<U’V’A)  = 0 



lv  fM  + 11273  42/3.1  < (1  + 1 ->2/3  ..2/3. 

7V-N,j+ltC=>  + 1)  -3  )+7V-N+l,j((j+1)  -<3-D  > 


+ V-N,j((^+1)2/3+^2/3"2C^~1)2/3)+v-N+l,j-l(32/3"(j-1)2/3) 


(4b) 


+ kN,i-i(i2/3-(i-1)2/3) 


Kt( A-N+l ,j“A-N, j )((3+l)2/3-(j-l)2/3) 


+ (A-N+i,j-i-A-N,:-i)(j2/3-(^1>2/3> 


3A 


N , j 


J^(u,v,A)  = 0 


,3)2/3 

V4 


T73  f (v-N+1 


2/3  . . 2/ 3 . 


,rv-N,j)((j+1)  } 


(4c ) 


+ (v-N+i,j-rv-N,j-i)(j2/3'<j‘1)2/3)] 


• i 
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b)  mesh  points  on  the  parabolic  line:  j = 0 


i = -N+l, . . . ,N-1 


Ui,0 


and 


9 J 


9 A 


v.  are  given 

1 5 U 


(u  , v , A ) = 0 =* 


1*0 


2.(|)2/3 

7T7T~  C(vi,l"vi-i,l)+(vi+i,o-vi-l,0>3 


= 3(u.  -,+u.  , , ) - ( 4u . n + U.  , n+U.,,  n) 
1,1  1-1,1  1,0  1-1,0  i+l,0 


This  scheme  is  consistent  and  the  order  of  the  equation  is 


Uy  * % * °«>2/3> 


i 

f 
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5 . Cauchy  problem  for  the  Tricomi  equation-Local  Formulation 

A ) Triangular  Isoparametric  Elements 

We  divide  the  domain  into  triangular  isoparametric 

elements  and  take  the  same  trial  functions  as  before. 

Let  E1’^  be  the  element  bounded  by  y = ^ and  the  two 

characteristics  and  S^’ ^ which  intersect  in  the  mesh  point 

(x  . • ,y . ) : 

i >3  3 


Step  j (j  = 1,...,2N) 

We  assume  that  the  solution  (u,v,A)  is  known  at  step  j-1, 
so  the  problem  is  well-posed  in  element  E1’-^  with  the  free  bound- 
aries  S^,J  and  Sj’.  We  take  the  variational  formulation  in 
every  single  element 


(u ,v , A)  = f [yu^-v^+A(u  -v  ) ]dxdy 

j J y x 

ri>  3 


A (udx+vdy ) 


i = -N , . . . , + N-J 


si»3us^  >: 


and  we  assume  that  A is  kept  fixed  on  the  free  boundaries. 


, ■ 
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3u . 


-J1,J  (u,v,A) 


= 0 - 


,3.  .4/3  1 8 r ,.7/3 
(iTh)  •— [ui,j(3  ' 


3 ,.  . . 10/3  6-10/3  27  ,.  n . 13/3 .27 . 13/ 3 

5^"i;  ~P  Tr3-1'  +rr3 


+ (ui,o-i+ui+i,:-i)(I(j-1)7/3+f(5-1)10/3+Tijl°/3+OT^-1)13/3 


27 .13/3 
130: 


)] 


Cla) 


J1,3(u,v,A) 

Ld 


o - 


vi}.(2j5/34(j-l)8/3-3j8/3-f|(j-l)11/3+fZj11/3) 

•Cvi,j-l  + Vi  + l,j-l)((^1)5/3+f(^1)8/3+^8/3  + ^(^1)11/3-^11/^ 


, 5/3 .. 5/ 3 3 .8/3.3, 


8/3. 


(lb) 


4”; — J^’3(u,v,A)  = 0 -* 

hi 


u . . , +u . 


(3)2/3 


u.  . . jjj-1  i-^-d-1  = 18,V  (v  _v  )(i5/3+(i-i)5/ 

i,3  2 5 hl/3  ^Vi  + 1 , j-1  i,j-lM;i  ^ 1; 

♦ f(j-l)8/3  - ^i8/3)  do) 

The  solution  of  step  j-1  gives  immediately  ^ (lb)  and  j (lc) 

and  substituting  u . into  (la)  we  get  A.  .. 

i » J 1 » 3 
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Step  0 : the  parabolic  line 


u.  „ and  v.  „ are  given, 
x , 0 i,0  & 


A.  n is  not  known. 
1 , u 


-gy- J1,1(U,V,A)  = 0 


i = - N , . . . , N - 1 


M 


B ^ 


(l,o) 


U - Ui  , 0+Ui  + l , 0 _ 27_  (4h)2/3  ( 

i,l  2 “20  h 1 / 3 Vi+1 , 0-Vi ,0 


On  the  other  side  (lc)  for  j = 1 will  give 


9Ai,l 


J^’^(u,v,A)  = 0 


Ui , 0+ui+l ,0  9 (fh)2/3 

2 10  “1/3 

h 


( v . . , „-v . n ) 

l + l ,0  i,0 


(lii) 


Ui  0 Ui  + 1 0’  v i 0’  Vi  + 1 0 are  *cnown»  so  and  (lii)  give  a 

different  value  for  u.  , and  A . n is  still  unknown.  So  we 

i,l  1,0 

cannot  take  the  variation  of  A on  the  parabolic  line  and  must 

get  the  values  of  A.  _ in  another  way: 

l , u 

The  analytic  connection  between  A arid  v on  the  parabolic  line 


2V  = Ax. 


So  we  can  get  the  values  o'  A.  _ by  numerical  or  an ‘lytic  inte- 

1 » V 

gration,  because  we  know  the  values  of  v.  n. 
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B ) Linear  elements 

For  j * jg  the  elements  EJ ,J  will  be  as  before  (triangu- 
lar isoparametric  and  the  trial  function  linear  in  £j,n- 

For  j £ the  elements  EJ’-'  will  be  ordinary  triangles 

bounded  by  y = y^_^  and  the  linear  approximation  of  the  charac- 
teristics and  Sj’-5  which  intersect  in  the  mesh  point 

( x . . ,y . ) . 

1 


Step  j (j  = 1»...»3q) 

We  assume  that  the  solution  (u,v,A)  is  known  at  step  j-1. 
We  take  the  variational  formulation  in  every  single  element 


J1, j(u,v,A)  = | [yu^-v^+ A (Uy-v^ ) ]dxdy  + j A(udx+v 


dy ) 


gi.luSj’3 


i = -N  , . . . ,+N-j 
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and  we  assume  that  A is  kept  fixed  on  the  free  boundaries, 


3u  . 


J 


(u , v , A ) = 0 


/ 3,  , 2 / 3 ,.2/3  ,,2/3,  lr  ,c.2 /3ih,.  ,,2/3, 

C^h)  •(]  -(]-!)  ) *IOtui,  j (6:  +1+(D-1)  ) 


+ ^i,j-l+ui  + l,j-lH2j2/3  + 3(j-l)2/3)] 


A . . , +A . , . , 

= . - ^J--1 3-IA.iJ-.1 

2 


(2a) 


9 v . . 
Li2 


■ J 1 ’ -1  ( u , v , A ) = 0 


2v.  . + v . . , + v . , 

1,3  l»3-1  i+ 1 j 3 - 1 


hai+l,j-l"Xi,3-l) 


(2b) 


yy  J1,J(u,v,A)  = 0 
"..-10 


Ui  , )-ltUi  + l ,9-1 


,3,2/3 

'V  ,.2/3 


2/3. 


1 5 3 


T7J  ‘1  )(vi+l,j-l-vi,j-l) 

(2c) 


We  get  again  3 explicit  equations  for  u.  v.  A.  .. 

i > 3 i>3  1 > 3 
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Step  0:  The  parabolic  line 


g and  g are  given 


i = -N, . . . ,N-1 


3 A 


3 - i 1 

J ’ ( u , v , A ) = 0 ■* 


LlL 


u . 


i >1 


ui,0»ul+1.0  = (!>2/3 


h 


T7T 


(v.,,  n_v-  n ) 

1+1,0  1,0 


On  the  other  side  (2c)  for  j = 1 will  give: 


3X7 


7 + » 1 


(u,v, A ) = 0 


i , 1 


u _ Ui,Q  + Ui  + l,0  _ iA)2/3 

i,l  2 " , 


(v.^-,  „-v,  n) 


173  vvi+l,0  vi,0 


( 2 i ) 


( 2 i i ) 


which  is  the  same  scheme  as  (2i). 

yy; J1’1  and  yy- — J1’1  give  the  same  scheme  because  the  trial 

i , 0 i,l 

functions  are  linear  in  x,y  in  every  element  so  u^  and  v^  are 
constant  in  every  triangle. 

In  order  to  get  the  values  of  A.  n we  integrate  the  equation  2v=A 
as  in  the  case  of  the  isoparametric  elements. 
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Numerical  Results 

In  the  step-by-step  method  we  must  store  only  3*(2N+1)  values, 
where  2N+1  is  the  number  of  points  at  the  parabolic  line.  So 
we  can  use  very  small  values  of  h and  get  any  accuracy  of  the 
solution  we  want  in  condition  that  the  solution  converges. 
Numerical  examples  are  given  in  figures  2,  3 where  the  results  of 
the  following  methods  are  compared: 

1)  The  analytic  solution 

00  , 3n+l 

Figure  2:  u(x,y)  = sinhx(y  + [ ? 3n+'l ) 3n<3n-2  ).  . ■ 3 j 

n=  1 

00  3n 

v(x,y)  = ccshxtl  * £ ifin-;).....!1 


y 3n 

(x,y)  = 2sinhx(  1 + [ fnTTn^~~  4 . 3 j 


Figure  u(x,y)  = 6x 

o 

v(x,y)  = 3y 
A(x,y)  = exy" 


2)  Solution  of  the  Cauchy  problem  in  the  total  hyperbolic  region 
where  the  elements  are  isoparametric  rectangles  bounded  by  the 
curves  £=const  and  n=const  and  one  row  of  isoparametric  triang- 
les adjacent  to  the  parabolic  line,  and  the  trial  functions  are 
bilinear  in  £ ,n  in  the  isoparametric  rectangles  and  linear  in 
£,n  in  the  isoparametric  triangles.  This  method  is  demonstrated 
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by  Yani v [ 1 ] . 

1 Step-by-step  method  with  isoparametric  triangles  and  trial 
functions  which  are  linear  in  £,n  i.e.  = 0. 

4)  Step-by-step  method  with  linear  corrections  near  the  parabolic 
line.  j q - 1 was  taken,  i.e.  in  one  row  the  trial  functions  are 
Linear  in  x,y  and  in*  the  others  linear  in  £,n* 

5)  Step-by-step  method  with  trial  functions  linear  in  x,y  in  the 
whole  hyperbolic  region,  i.e.  jg  = 2N. 

Initial  values  are  given  for  u,v  and  A in  all  the  methods  des- 
cribed above. 

In  both  examples  every  second  step  is  noted  ( j-Q  ,2  ,4  , . . . ,12) . 


V = 2y 
A = 6 xy 2 

The  table  shows  the  numerical  results  for  v for 
different  values  of  h in  the  step-by-step  method 
with  trial  functions  which  are  linear  in  x and  y, 

Convergence  0(h)  is  demonstrated 
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2 

v-analytic=  3y 

Cauchy  in 
every  triangle 
linear  in  £,n 

Cauchy  in 
every  triangle 
linear  in 

K >n(x,y  for  j=l) 

Cauchy 

every 

angle 

ear  in 

0 

0 . 000 

0.000 

0 . 000 

0 .000 

1 1 

0 . 1875 

0 .000 

0 . 000 

0.000 

2 

0.472 

0.395 

0 .384 

0.375 

3 

0.811 

0.564 

0 . 565 

0.570 

4 

1.191 

1.075 

1.064 

1.053 

5 

1.603 

1.320 

1.321 

1.329 

6 

2 . 044 

1.902 

1.891 

1.878 

7 

2.511 

2 . 200 

2.201 

2.211 

8 

3.000 

2 .837 

2.825 

2.811 

9 

3.510 

3.178 

3.178 

3 . 190 

10 

4. 040 

3.858 

3 . 846 

3 .831 

11 

4 .587 

4.235 

4 . 236 

4 . 248 

12 

5.151 

4.954 

4.942 

4 .926 

Table  .2:  Numerical  results  of  v for  j = 0 , 1 , 2 , . . . ,12  , 

h = k 

analytic  solution:  u = 6x 

o 2 

v = 3y 


A = 6xy‘ 


Nonuniform  (i.e.  oscillating)  convergence  is  seen. 
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FINAL  REMARKS 

Tricomi's  problem,  with  a view  of  application  to  transonic 
flows  and  nonlinear  generalizations,  has  induced  much  ingenious 
effort,  with  the  progress  attained  ever  increasing.  A few  prob- 
lems, however,  need  further  clarification.  Unresolved  problems 
exist  regarding  the  local  behavior  of  approximations  near  the  para- 
bolic line,  the  local  behavior  of  weak  and  strong  discontinuities 
(is  a local  loss  of  accuracy  inherently  unavoidable?)  and  the  pos- 
sible pollution  of  the  whole  field  by  local  errors  at  boundaries 
and  non-smooth  regions.  Further  clarification  is  needed  with  regard 
to  observed  local  numerical  phenomena  (e.g.  decoupling  of  points, 
local  under  or  super  convergence  and  accuracy)  and  to  a better 
understanding  of  numerical  procedures  at  the  parabolic  (sonic) 
and  discontinuity  (shock  waves)  lines.  Possible  improvements  are 
sought,  trying  alternatives  to  finite  difference  methods  (e.g. 
Murman-Cole)  and  classical  weak  solutions  (e.g.  Friedrichs)  via 
new  formulations  as  the  one  proposed  here,  and  finite  elements 
following  the  inherent  structure  of  themathematical  field  (grid 
of  characteristics,  where  real). 

The  results  above  are  rough  and  exploratory.  In  view  of 
new  ways,  further  work  and  the  generality  of  the  underlying  formu- 
lations, both  analytic  (variational)  and  discrete  (characteristic 
grid  with  local  treatment  where  needed)  are  indicated,  addressing 
the  questions  above.  Extensions  at  work  are  to  3-dimensional  lin- 
ear Tricomi-type  problems,  and  to  nonlinear  mixed  type  equations. 
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1 . Introduction 


Systems  of  nonlinear  conservation  laws,  of  the  form 


(1.1) 


u.  + f(u)  = 0 

t X 


with  u,f  n-dimensional  vectors,  are  said  to  be  hyperbolic 
in  a region  H (of  u-space)  if  for  every  u£H,  the  eigen- 
values of  fu(u)  are  real  ami  distinct.  For  example,  in 
the  system 


(1.2) 


p.  - q =0 
Hx 


;qt  - (p  )x  = 0 , 


H is  the  half-plane  p > 0.  The  system  (1.2)  will  be  adop- 
ted on  a'  prototype  problem  for  the  present  discussion,  as  it 
arises  in  various  applications  [14,  15]..  In  the  case  p > 0, 
the  uniqueness  of  the  solution  of  the  Cauchy  problem  for 
(1.2)  is  obtained  in  [15]. 

Some  systems  of  conservation  laws  admit  an  entropy  func- 
tion, i.e.  a convex  scalar  function  U(u)  such  that  [ 6 , 8 J 


(1.3) 


f XU  = F 
u u u 
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for  another  scalar  function  F.  This  is  true  in  particular 
for  pairs  of  equations  and  for  some  larger  systems,  such  as 
the  equations  of  compressible  fluid  flow.  Multiplying  (1.1) 
by  Uu  and  using  (1.8),  it  follows  that  smooth  solutions  of 
(1.1)  satisfy  an  additional  conservation  law 


(1.4) 


U ( u ) . + F( u ) = 0. 

L X 


For  a given  system  one  typically  looks  for  'functions  U,F 

satisfying  (1.4).  Alternatively,  a differential  equation  for 

U may  be  derived  from  (1.3).  For  example,  multiplying  the 

2 

first  equation  in  (1.2)  by  p and  the  second  by  q,  it 
follows  that  (1.4)  holds  with 

3 2 9 

(1.5)  U( p ,q ) = P /3  + ^ /2,  F(p,q)  = -p  q. 


We  emphasize  that  (1.4)  does  not  hold  in  general  for 
weak  solution  of  (1.1).  Nevertheless,  there  are  (for  our 
purposes)  two  important  implications  of  the  existence  of  an 
entropy  function.  One  is  the  existence  of  a change  of  vari- 
ables, given  by  [2] 


(1.6) 


v = Uu(u) 
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which  puts  the  system  (1.1)  into  generalized  symmetric  form, 
[3]: 


(1.7) 


(4>v(v))t  + OP  (v))  = 0 


where  4’,’!'  are  scalar  functions  of  v.  In  fact,  4> , V are 
given  explicitly  by 


(1.8)  <P(v)  = u-Uu(u)  - U(u ) , f(v)  = f(u)*U  (u)  - F(u); 


assuming  that  all  the  indicated  derivatives  exist,  then  diff- 
erentiation of  (1.8)  shows  that 


(1.9) 


u = <t>v(v)  , f (u)  = (v) 


which  establishes  the  equivalence  of  (1.1)  and  (1.7),  even 
for  weak  solutions.  For  our  example  (1.2,  1.5),  we  readily 
obtain 

r\j  O.  T 2 T 

(1.10)  v = (p,q;  = (p  ,q)  , 


(1.11) 


, 2%3/2 

4>(v)  = jp 


]>2 

2^ 


Y(v) 


-pq> 


and  from  (1.7),  the  corresponding  symmetric  system  is 
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(1.12) 


'X. 


0 


0 . 


Even  for  smooth  initial  data,  the  solution  of  (1.1) 
is  not  unique;  the  convexity  of  U is  used  as  a means  of  re- 
moving the  nonuniqueness.  We  suppose  that  u,  a weak  solu- 
tion of  (1.1),  is  the  limit  as  e -*  0+  , boundedly  in  meas- 
ure (two-dimensional  measure  in  the  x,t  plane)  of  the  solu- 
tion u of 
e 


(1.13) 


ut  + f(u)x  = euxx, 


the  solution  of  which  is  smooth.  Multiply  (1.13)  by  $U  , 

00  , 
where  s€Cq  is  a nonnegative  test  function  of  x,t,  and 

integrate  over  x,t.  After  partial  integration  we  obtain 


U(ueHt  + f(ue)4)x 


= -e 


e ,x 


(U  (u  )u  H + e 


uu 


e ,x 


U ( u ) 4>  . 

e xx 


Passing  to  the  limit,  we  obtain,  using  U > 0, 


(1.14) 


U(u)t  + F(u)x  < 0 


in  the  sense  of  distributions.  When  f is  genuinely  non- 
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linear  [7]  the  entropy  inequality  (1.14)  is  equivalent  to 

the  classical  entropy  condition  for  weak  solutions  of  (1.1) 

[9],  at  least  in  the  small;  the  theory  is  extended  to  include 

strong  shocks  in  [11,  13].  Such  entropy  inequalities  are 

readily  incorporated  into  difference  schemes,  including  such 

schemes  to  automatically  select  the  physically  correct  weak 

solution  [5,  8,  10,  12,  13].  ^ 

The  convexity  of  U is  also  related  to  the  hyperbolic-  ; 

ity  of  (1.1),  as  shown  by  the  following  theorem.  j 

Theorem  1.1:  Suppose  that  the  system  (1.1)  admits  an 

entropy  function  U,  convex  for  u€H.  Then  the  eigenvalues 

of  fu(u)  are  real  for  all  u€H. 

Proof . Differentiating  (1.6)  with  respect  to  u and 

(1.9)  with  respect  to  v,  we  have  U ^(u)  = $ (v),  where 

uu  vv  ’ 

v = U (u),  so  that  the  local  convexity  of  U with  respect 
to  u is  equivalent  to  that  of  <I>  with  respect  to  v.  The 
eigenvalues  A of  f^(u)  are  also  the  eigenvalues  of  the 
mixed  symmetric  problem. 

(1.15)  <t>  r = A4>  r 

Vv  vv 

where  r is  the  corresponding  eigenvector.  Since  V is 
symmetric  and  4>vv  symmetric  and  positive  definite,  the  eig- 
envalues A are  real.  Thus  when  hyperbolicity  fails  in  (1.1), 


i 
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ior  example  we  allow  p < 0 in  (1.2),  the  convexity  of  U 
fails,  and  (1.14)  no  longer  follows  from  (1.12).  In  fact, 
we  may  question  whether  (1.12)  is  an  appropriate  form  of 
regularization  in  this  case.  Even  in  the  purely  hyperbolic 
case,  it  is  known  that  different  choices  of  regularization 
can  lead  to  different  weak  solutions  as  limits  [1,  4,  16]. 

We  can  recover  (1.14),  even  with  U not  necessarily  convex, 
by  using  a regularized  form  of  (1.7).  Indeed,  let  v be 
the  limit  of  v^ , satisfying 


(1.16) 


<w>t 


(4'  (v  )) 
v e x 


= e(Avr.  ) 
e ,x  x 


where  A is  positive  definite,  symmetric;  A may  even  de- 
pend on  v£ . Multiplying  (1.16)  by  v&  and  integrating, 
we  use  (1.6,  1.8,  1.9)  to  obtain 


0.17) 


U(u)  = v*<&v  - 4> , F(u)  = v • T v~4' 


and  readily  obtain  (1.14)  as  above.  However,  when  hyper- 
bolicity  fails,  the  transformation  u-*-v , described  in  our 
sample  problem  by  (1.10),  is  in  general  no  longer  invertible, 
so  that  (1.1)  and  (1.7)  are  no  longer  equivalent. 

These  problems  may  be  aggravated  in  general  by  the  existence 
of  many  entropy  functions,  each  with  its  corresponding  gen- 


— * 


f 

I 

i 
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eralized  symmetric  system. 

In  order'  to  approximate  the  physically  relevant  solutions 
of  a given  system  (1.1)  of  mixed  type,  we  require  knowledge 
of  how  much  accitional  physics  is  needed  to  make  the  problem 
well-posed,  in  addition  to  the  specification  of  proper  initi- 
al/bounda  y conditions.  For  example,  which  dependent  varia- 
bles should  be  used,  which  form(s)  of  regularization  are  app- 
ropriate, and  are  entropy  inequalities  (1.14)  sufficient? 

1 ! 

Obviously,  these  questions  are  not  independent. 

In  section  II,  we  briefly  review  the  theory  of  piecewise 
smooth  weak  solutions,  in  slightly  generalized  form.  In  Sec- 
tion III,  we  apply  these  methods  to  a class  of  problems  inclu- 
ding (1.2),  and  obtain  suitable  forms  of  regularization  and 
entropy  conditions.  These  results  are  incorporated  into  a 
difference  scheme  for  (1.2)  in  Section  IV. 

i I 


I 
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II.  Discontinuous  solutions  of  systems  of  mixed  type 
We  next  consider  systems  of  the  form 


(2.1)  g ( u ) + f ( u ) = 0,  u€D  = HUEUOHnBE) 

x y 

which  are  hyperbolic  in  H only.  For  example,  we  rewrite  our 
model  problem  (1.2)  in  the  form 


(2.2) 


(P2)x-qy  = 0 
qx-py  = 0. 


Indeed,  we  have  in  mind  primarily  the  case  that  (2.1) 
is  a pair  of  equations,  and  is  elliptic  in  E.  Thus  in 

(2.2),  H(E)  is  the  half-plane  p > 0 (p  < 0). 

A system  (2.1)  will  admit  an  entropy  function  if  there 
exists  a vector  w(u)  such  that 

(2.3)  g^(u)w(u)  = U (u),  and  f^(u)w(u)  = F (u). 

u u ’ u u 


As  previously,  the  existence  of  an  entropy  function  im- 
plies that  the  system  can  be  put  in  generalized  symmetric 
form;  setting 


50 


(2.4)  <Mw)  = g(u)  •w(u)-LKu)  , t(w)  = f(u)*w(u)  - F(u), 

it  follows  from  (2.3)  that 

(2.5)  4>w(w)  = g(u)  and  4*w(w)  = f(u), 
so  that  (2.1)  is  equivalent  to 

(2.6)  (*  (w))  + (¥,  ,(w))„  = 0 

w x w y 

as  long  as  the  mapping  of  D+R  determined  by  u-*w(u) 
invertible.  The  convexity  of  U (with  respect  to  u) 
of  $ (with  respect  to  w)  are  no  longer  equivalent, 
ever. 

Our  example  (2.2)  is  already  in  symmetric  form, 

(2.7)  u = (p,q)T  , 4>(u)  = p3/g  + q2/2  > 'P(u)  = -pq. 

Thus  we  may  set  u = w and  find  U,F  from  (2.4 
taining 


(2.8) 


U(u)  = jp3  + ^q2  , F(u)  = -pq. 


is 

and 

how- 

with 


) , ob- 


However,  there  are  other  possible  choices,  for  example 
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(2.9)  w = (p3+|q2 


„ 2 . T 
3p  q)  , 


U(u)  = |-p5  + |-p2q2,  F(u ) = 


3 13 

-p  q-^i 


Corresponding  to  this  choice  of  w,  there  is  a symmetric 
form  (2.6),  with  <t>(w)  , 'J'(w)  obtained  from  (2.4,  2.9). 

Thus  for  the  pair  (1.2)  we  have  three  possible  entropy  func- 
tions: fr,  given  by  (1.5);  U,  given  by  (2.8);  and  U,  given 
by  (2.9).  More  could  be  constructed.  U and  ll  are  convex 
in  H;  U is  convex  when  p * (^q  ) . To  each  entropy  fun- 

cion  corresponds  a choice  of  dependent  variable  (1.10),  (2.7) 
of  (2.9),  and  a corresponding  symmetric  system,  two  examples 

of  which  are  (1.12)  and  (2.2).  These  systems  are  all  equi- 

3 2 1/3 

valent  in  the  region  p > ( jj-q  ) ; the  two  systems  (1.12) 

and  (2.2)  are  equivalent  for  u£H,  but  not  for  all  u€D. 

To  decide  which,  if  any,  of  these  systems  can  correspond  to  a 
well-posed  problem,  we  consider  weak  solutions  containing 
shocks . 

For  a specific  choice  of  dependent  variable,  denoted  by 
w,  we  consider  our  system  in  symmetric  form  (2.6).  Let  us 
assume  that  w has  an  isolated  discontinuity,  along  a shock 
curve  in  the  x,y  plane,  described  by 


(2.10)  -sinGdx  + cosOdy  = 0. 


Equation  (2.10)  corresponds  to  s = tanO , where  s is 


IV 


liL 


I 
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the  shock  speed,  in  the  case  where  x is  a time  variable. 
However,  in  systems  of  mixed  type,  infinite  shock  speeds  are 
possible;  thus  the  polar  representation. 

The  Rankine-Hugoniot  relations,  expressing  conservation 
of  material,  are 

(2.11)  -sin9[$>  ] + cosStT  ] = 0 s 

w w 

where  [ ] denotes  the  jump  in  a quantity  across  a shock,  as 

usual.  At  any  point  w,  the  characteristic  speeds  tan 

and  corresponding  eigenvectors  r^ , i = l,...,n,  are  given  by 

(2.12)  cosn  • f r.  = sinn-  $ r.  ; 

the  system  (2.6)  is  hyperbolic  at  w if  the  are  real 

and  distinct,  and  elliptic  at  w if  the  are  all  com- 

plex. The  system  (2.6)  is  hyperbolic  whenever  $ is 
locally  convex.  Assuming  that  3>  is  locally  convex  for  all 
w€H,  and  that  H is  a convex  region  (of  w-space),  shocks  of 
infinite  speed,  corresponding  to  cos0  = 0,  [<t>w]  = 0,  cannot 
connect  two  points  in  H,  but  such  shocks  might  connect  a 
point  in  H with  one  in  E. 

Given  a point  w^ , let  r(wg)  denote  the  set  of  points 
in  D which  can  be  connected  to  w^  by  a shock,  i.e.  for 
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which  the 
some  value 
wu  itself 
de  ci i bed 


Rankine-Hugoniot  relations  can  be  satisfied,  for 
ol  0.  By  definition,  Kw^)  does  not  contain 
. I'(wq)  is  a family  of  one-parameter  manifolds 
by  the  differential  equation 


5 


(2.13)  cosoCcosoH1  (w)-sino$>  (w))w’  = o ’ (4>,.(w) -4>„(wn  ) ) , 

WW  Ww  w w 0 ’ 


where  primes  denote  differentiation  along  T(wq),  and 
o = o(w,wQ)  is  the  value  of  0 for  which  the  Rankine-Hugon- 
iot relation  (2.11)  holds  between  w^  and  w.  We  next  ass- 
ume that  [11]. 


(2.14) 


o'  / 0 for  all  Wq€D,  w6T(wq) 


where  o'  is  determined  from  (2.13),  with  w'  > 0.  This 

is  equivalent  to  the  assumption  that  tana(w,WQ)  is  never  one 
ol  the  characteristic  speeds  at  Wq  or  at  w.  For  w€E, 
this  condition  is  always  satisfied,  as  the  characteristic 
speeds  are  complex.  Thus  (2.14)  is  really  only  an  assump- 
tion about  the  hyperbolic  region. 

The  condition  (3.14)  was  postulated  in  [11]  as  a 
statement  of  strong  nonlinearity  of  the  system  (1.1),  assumed 
to  be  hyperbolic.  When  (1.1)  or  (2.1)  is  hyperbolic, 
this  condition  implies  the  classical  condition  of  genuine  non- 


I 

I 


I 
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linearity  1 7 J , and  is  almost  equivalent  for  pairs  of  equat- 
ions. For  Wg6H,  r(wQ)nH  consists  of  2n  smooth  distince 
one-parameter  manifolds,  each  with  Wg  as  one  of  its  limit 
points,  and  the  other  limit  point  on  the  boundary  of  H or 
at  infinity.  Near  wQ , the  manifolds  of  r(Wg)  are  in  the 
direction  of'  ±r.(Wg),  depending  on  whether  o increases  or 
decreases  as  one  moves  away  from  Wg. 

For  Wg€E,  it  follows  from  (2.13)  that  r(Wg)  does 
not  contain  Wq  as  a limit  point,  for  if  it  did,  at  least 
one  of  the  characteristic  speeds  at  Wq  would  have  to  be  real. 
Thus  there  are  no  sufficiently  weak  shocks  connecting  two 
points  in  E. 

For  the  system  (2.1),  the  jump  conditions  (2.11)  are 


(2.15) 


where  w 
obtain , 


2 2 

sin0(p  ~Pq)  + cos0(q-qg)  = 0 

sin0(q-qg)  + cosO(p-pg)  = 0 
T T 

(p,q)  , wQ  = (pg,qg)  . From  (2.15),  we  readily 


I 


I 

( 


(2.16) 


2 2 
sin  0(p+pg)  = cos  0, 


so  that  p+ Pg  i 0,  i.e.  there  are  no  shocks  connecting  two 


points  in  L. 


in  cliff  erent 


reg  1 on; ; . 


Below  we  show  I’(Wq)  for  w[i 


I 


T(Wq)  for  wQ6E  (p0<0) 


Corresponding  to  the  symmetric  system  (2.6),  there  exist! 


an  entropy  function 


(2  I /) 


!J(w)  = w<t  — <t> , F(w)  = w • *P  - 4* : 

w w 


the  entropy  inequality  (1.14)  is  equivalent  to  regaining 
that  ( 8 ] 


(2.18) 


-sinO[U]  + cosO[F]  £ 0. 


For  shocks  which  are  limits  of  viscous  profiles,  obtain- 
ed by  regularizing  (2.6)  as  in  (1.13)  or  (1.16),  (2.18) 
holds  with  strict  inequality.  In  the  simplest  case,  when  the 
right  side  of  (2.6)  is  replaced  by  a second  directional  de- 
rivative, the  viscous  profile  between  two  states  w , w+ , is 
the  solution  of  the  system  of  ordinary  differential  equations 


- -sino  (w.  ,wn  ) (<t>  (w(z)  ) -0  (wn  ) ) + 
dz  1 0 w w 0 


(2.19) 


+ coso(w.,w  )(f  (w(z))-4'  (w,)), 

T — W W T 


W ( +°°  ) = W+  , W ( — 00  ) 


Let  us  examine  the  implications  of  (2.18)  as  applied 
to  the  system  (2.2),  utilizing  each  of  our  entropy  functions. 
With  U,F  given  by  (2.8),  we  use  the  jump  conditions  (2.15) 

to  obtain 
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( 2 - ? 0 ) -sin0[U]  + cos0[F]  = -|-sin0  (p-pg ) 3 

which  is  the  entropy  jump  in  a shock  connecting  two  states 
(pgj'lg)  and  (p,q).  Thus  for  this  entropy  function,  we  see 
that  (2.18)  implies  that  either  that  portion  of  r(pg,qg) 
to  the  right  or  to  the  left  of  the  line  p = Pq  is  admissible, 
depending  on  whether  (p^q^)  is  the  upstream  or  downstream 
point . 

The  entropy  function  (1.5)  was  obtained  with  the  vari- 
ables x,y  switched,  which  is  equivalent  to  switching  ft  and 
Therefore  to  compute  the  entropy  jump  in  this  case,  we 
compute 

(2.21)  -sin0  Lh  + COS0  [&]  = I sin0(p-po)2(q-qg). 

In  this  case,  either  that  portion  of  r(pg,qg)  above 
or  below  the  line  q = qQ  is  admissible,  and  possibly  the 
shock  corresponding  to  p = “Pg » q = qg  . 

Finally  for  the  entropy  function  L? , given  by  (2.9),  we 
obtain  (in  the  special  case  qg  = 0) 

(2.22)  -sin0[U]  + cos0[F]  = -ig-sinO  (p-pQ  ) 3 (4p2  + 2ppQ-pg  ) , 

and  the  curve  r.(p  ,qg)  is  divided  as  shown  in  Fig.  2,  the 
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III. 


(3.1) 


(3.2) 


(3.3) 


(3.4) 


A class  of  pairs  of  mixed  type. 

We  consider  the  class  of  pairs  of  equations 

U'  (p)x-qy  = o 


T 

which  corresponds  to  (2.6)  with  u = (p,q) 

2 

1>(p,q)  = Up)  + q /2>  Up,q)  = -pq 

and  admits  a corresponding  entropy  function 

2 

U(p,q)  = u*1>u-<t>  = pU(p)-Up)+q  /2>  F(p,q)  = -pq. 

As  in  the  sample  problem  (1.2),  which  is  of  course  a 
special  case  of  (3.1),  there  exist  other  functions  ft,  ? 
such  that  (1.14)  holds,  for  example 

2 

U(p,q)  = qU(p),  F(p,q)  = -Up)  - q /2 
The  system  (3.1)  is  hyperbolic  whenever  £"(p)  > 0; 


the  specific  assumptions  we  make  on  C(p)  are 


I 


-60- 

5"(p)  > 0,  p > 0,  5" (p)  < 0,  p < 0, 

V " (p)  > 0,  p > o, 

the  assumption  (3.6)  assuring  that  the  system  is  strongly 

nonlinear  in  the  sense  of  [11]  in  the  hyperbolic  region. 

Without  loss  of  generality,  we  take  5(0)  = 5'(0)  = 0. 

The  Rankine-Hugoniot  relations  (2.11),  describing  which 

states  (p,q)  can  be  connected  to  a given  state  (pg,qg) 

are  in  this  case 

I 

( 

(3.7)  sin0(5'(p)  - 5'(pg))  + cos0(q-qQ)  = 0. 

(3.8)  sin0(q-qQ)  + cos0(p-pg)  = 0. 

From  (3.7,  3.8),  it  is  clear  that  there  are  no  shocks 
with  0=0,  but  shocks  between  the  elliptic  and  hyperbolic 
regions  with  Q - ± may  exist.  From  (3.2,  3.8),  we  easily 

obtain 

5'(P)  - 5 ' (p0 ) 2 

(3.9)  — = cot  0 * 0; 

P-P0 

since  the  elliptic  region  (p<0)  is  convex  and  5,<:0  there, 
we  have 


(3.5) 

(3.6) 


i 


61 


Theorem  3.1:  There  are  no  shocks  connecting  two  states 
in  the  elliptic  region. 

Let  us  therefore  take  (p0,qQ)eH,  qQ  = 0 without  loss 
,,  generality,  and  ask  about  the  shape  of  the  curve  T ( pQ  ,qQ ) . 
from  (3.6),  we  know  that  6'  i 0 (prime  denoting  differen- 
tiation along  D on  each  branch  of  r(pg,qg).  Then  from 
(3,8)  we  have 

Theorem  3.2:  Each  branch  of  of  Kp^q^)  is  star- 
shiped  about  the  point  pQ,q^. 

It  follows  from  (3.6,  3.8,  3.9)  that  for  p < pQ 

i 

I 

q-cin 

(3  10)  cote  = y.  < (^"(p-))1^ 

1 P"Po  0 

and  then  from  Theorem  3.2  and  (3.10) 

Theorem  3.3:  For  (p^q^)  € H,  the  two  branches  of 
f(p0 ,qQ ) initially  in  the  direction  of  decreasing  p must 
enter  the  elliptic  region.  These  two  branches  will  eventual- 
ly join  (as  in  Fig.  1 above)  if  £'(p)  = ^'(pQ)  for  some 
p < 0.  If  E'(p)  < £'(pg)  for  all  p < 0,  the  two  branches 
will  remain  apart,  as  shown  in  Fig.  3 below. 
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Fig-  „.3.:  F(p0,q0),  (p0,qQ)  £ H,  5'(p)  < 5:(pQ)  for  all  p<0 

The  situation  described  in  Fig.  3 can  only  arise  for 
sufficiently  large  pQ ; for  small  positive  pQ  , r(p0,qQ) 
will  be  as  shown  in  Fig.  1.  For  (p0,qQ)€E  or  on  the 
boundary  between  E and  H,  the  curves  of  Kp^jq^)  res- 
emble those  of  Fig.  1. 

The  entropy  inequality  (2.18),  applied  to  (3.3),  (3.4) 
respectively,  gives 


(3.11)  -sin0[U]  + cos9[F]  = sin0[  £ (p ) -F,  (p^  ) - j(  p-Pg  ) ( £ ' ( P ) + £ ' ( PQ  ) ) ] , 

% a,  c<p)-ecpQ  > , 

(3.12)  -sin0[U]  + cos9[F]  = sin0(q-qn)[ - y(  £ ' ( p ) + £ ' ( pn  ) ) ] 


Exactly  as  in  the  sample  problem  (2.2),  the  entropy  func 
tion  U leads  to  a division  of  r(pg,qg)  with  respect  to 


I 

I 
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i 

p-pQ,  and  the  entropy  function  & to  a division  of  r(p0,qQ) 
with  respect  to  q— q q . Alternatively,  the  inequality  (1.14) 
as  applied  to  the  entropy  function  (3.3)  is  equivalent  for 
piecewise  smooth  solutions  to  the  inequality 


! 


i 

\ 


1 

1 

, 


t| 

►4 


(3.13) 


p < +00 ; 


( 1 . 14 )  with  the  entropy  function  (3.4)  is  equivalent  to 


(3.14) 


+00 . 


It  is  quite  natural  that  different  entropy  functions 
should  give  different  admissible  shocks,  as  they  correspond 
to  different  forms  of  regularization  of  (3.1).  The  system 
(3.1)  is  already  in  the  symmetric  form  corresponding  to  U, 
and  so  any  regularization  of  the  form  (1.16)  leads  to 

(1.14).  The  symmetric  form  of  (3.1)  corresponding  to  U 

• , ' r\j  % 

given  (3.4)  involves  the  transformation  p = £'(p),  q = q, 
which  is  not  invertible  when  p can  ^hange  sign.  Alternati- 
vely, we  could  obtain  (1.14)  for  & by  a regularization  of 
(3.1) 


(3.15) 


CCp’x-qy  = qxx 
qx-Py  « <C'Cp»xx 


! 


t 
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which  may  or  may  not  be  physically  reasonable.  Some  knowledge 
of  the  physical  dissipation  mechanism  is  necessary  to  choose 
the  correct  form  of  regularization  and/or  entropy  inequality. 
In  view  of  the  results  of  [1,  17]  it  seems  likely  that  this 
is  true  even  in  the  purely  hyperbolic  case.  We  conclude  this 
section  by  proving  that 

Theorem  3.4:  The  shocks  allowed  by  (2.18),  (3.3) 
(alternatively  (3.13))  are  the  limits  of  viscous  profiles. 

Proof : Let  (p+,q+),  (p_,q_)  be  two  states  which  can 

be  connected  by  a shock,  with  p+  > p_.  Let  9 be  such  that 
(3.7),  (3.8)  are  satisfied  (between  (p+>q+)  and  (p_,q_)), 
with  sin0  >0  for  definiteness.  We  wish  to  show  the  exist- 
ence of  a solution  of  (2.19),  which  becomes 


(3.16) 


Pz  = -sin0(£'(p)  - 5’(Pq))  - cos^(q-q+) 
qz  = -sin0(q-q+)  - cos0(p-p+) 


p(±»)  = p , q(±«)  = q 

+ + 

It  follows  from  (3.8)  that  q+-q  has  the  same  sign  on 
cos0,  which  we  take  positive.  Let  ft  be  the  rectangle  with 
sides  parallel  to  the  p,q  axes  and  (p_>q  )»  (p+>Q+)  at  oppo- 
site corners,  as  shown  in  Fig.  4.  From  (3.16)  the  orbit 


— * 
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leaves  ft  on  all  four  sides. 


Using  the  strong  nonlinearity  condition  (3.6),  it  is 
not  difficult  to  prove  that  (p+>q+)  is  an  improper  node  and 
(p  ,q  ) is  a saddle,  with  respect  to  the  vector  field  gen- 
erated by  (3.16)  [11].  This  holds  even  if  (p_,q_)  is  in 

the  elliptic  region.  It  is  also  easily  shown  that  there  are 
no  other  critical  points  in  ft,  so  that  the  orbit  entering 
the  saddle  from  within  ft  must  have  come  from  (p+,  q+)- 
This  completes  the  proof,  and  also  shows  the  uniqueness  of 
the  orbit  up  to  translation,  as  desired. 


I 
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IV.  A difference  scheme  for  the  small  disturbance  problem 

In  this  section  we  utilize  the  results  above  to  construct 
a finite  difference  approximation  to  (2.1),  which  automatic- 
ally excludes  unphysical  solutions  and  which  approximates 
(2.1)  to  second  order  in  the  mesh  size  even  when  the  solu- 
tion is  not  smooth.  The  geometry,  boundary  conditions,  and 
assumed  form  of  the  solution  are  shown  in  Fig.  5, 


. 


in  which  Pq  > 0 and  the  smooth  function  q^(x)  are  given. 
We  assume  that  the  solution  is  piecewise  smooth,  and  of  boun- 
ded variation  in  x,  uniformly  in  y.  The  shock  curve  will 
reach  the  x-axis  in  general,  but  we  assume  that  it  does  not 


reach  the  lines  x = L or  x = 0. 


The  discrete  variables 
shown  in  Fig.  6.  The  q's 


k k+1/2  . . . 

p.,  q . are  oriented  as 

1 H]+l/2 

are  oriented  at  the  0-points 


I 


I 


I 


is  our 


and 


q 


j + 1/2 


approximation  to  q(x,y),  for  x,y  in  the  rectangle 
j Ax<x< ( j + 1 ) Ax , kAy<y< (k+1 )Ay . The  p's  are  oriented  at  the 
x-points  in  Fig.  6;  however,  their  interpretation  is  more 
complicated.  Let 


5 + 1/2 


1/  k k . 

7<pj+Pj+i>  * 


then  in  the  rectangle  ( j - j)  Ax<x<  ( j +|o  Ax  , (k-i)Ay*y*  (k  + i-)  Ay  , 
our  approximation  to  p(x,y)  is  piecewise  linear,  given  by 


(j+i)Ax-x  . x-(j-y)Ax  , 

Ph(x>y}  = C h ]Pj-l/2  + C Ax  ] P j -1/  2 ’ 


. . k k+1 

The  difference  approximations  determining  the  p^ , Qj+2 


are  as  follows: 
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(4.3) 


(4  4 ) 


k + 1 / 2 _ k‘  1'2  k + 1 k 

— - -1— jrr1  = 0>  3 = 1,2,.  ..  ,M-1,  k = 0 , 1 , . . ,N-1  ; 


,k  2 , k.2  k + 1/2  k-1/  2 

iPi±l>  -(pj>  q,U/2-qj  + l/2 


= 0 , j - 0 , M- 1 , k • 0 , 1 , . . ,N; 


/ k k w k . k , k k+1/2  k-1/2 

( ( b ) Ml*3/2~Mi-l/2)(yi+3/2Mi+l/2yi-l/2  _ qi+l/2~qj  + 1/2 

3Ax  Ay 


lk  k 

'Ax^j  +5/2f'i+5/2“4Mj  + 3/2C  j + 3/2  + 


+ 6yj+l/2C j+l/2"4yj-l/2Cj-l/2  + Pj-3/2Cj-3/2  ] 


j = 1,2,...  ,M-2 , k = 0,1 ....  ,N. 

In  (4.5),  C j +1/2  = C ( ( j + 1/ 2Ax ) , £ a nonnegative  c” 

function  chosen  so  that  £ > 0 in  a region  containing  the  sh 
shock,  and  C = 0 near  the  boundaries  in  x,  so  that  (4.5) 
does  not  require  any  p-values  outside  the  rectangle. 

We  first  show  that  the  difference  scheme  (4. 3-4. 5)  sat- 
isfies a discrete  form  of  the  inequality  (1.14)  for  the 

]<- 

entropy  function  (2.8).  Let  <|>j  = <J>  ( j Ax  ,kAy ) , <J>  a nonnega- 
00 

tive  Cn  test  function;  multiply  (4.3)  by 


k+1/2  k+1/2 

a*  ay  (Vl/2^3-1/2)  ^1/2 
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k k 

and  (4.5)  by  Ax  Ay  ^ j + 1/2  ^j+l/2’  ac^  an<J  sum  over  j >k 


We  assume  that  the  discrete  variables  are  uniformly  bounded, 
and  use  the  smoothness  of  C and  the  test  function  <\>  re- 
peatedly. After  several  partial  summations,  we  obtain 


, ,k  ,3  . k .2  k . k , k ,2W  k .3 
AxAy  l 3 + 1/2  (Mi-H/2}  ^ j -1/2  M j + 1/ 2 ^ j -1/ 2 (llj-l/2) 

• 1 6 

] >k 


;*]>!/»•  S-l/2 


Ax 


, k+ 1/2  . k + 1/2 


1,  k+1/ 2 , 2 , i + 1 

* 2"tqj +1/ 2 J '-1 — sr — > 


k+1  k 

k+1/2, Mi+l/2  Mi+l/2, 
qj  + l/2(  2 } 


,k+l  *k 
^j  + l/2~^j+l/2  j 
Ay  ;J 


Ay  j^k(lJj  + 3/2(^i  + 3/2<<)j  + 3/2) 


1/2 


9,,k  (r  Ak  ,1/2.,  k ,r  . k ,1/2,2 

2pj  + l/2(cj  + l/2<J>j  + l/2)  yj-l/2(Cj-l/2^’j-l/2) 


+ 0 ( Ax+Ay  + Ax  sup  I | P j + 2 ~P j I ) 
k j 


The  left  side  of  (4.6)  is  an  approximation  of 


j ( (|p'3+|<l2Hx-pq<>'y)dxdy ; 


..  «w_.  *»  ■ . 
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the  right  side  is  a nonnegative  term  plus  a remainder  which 
can  be  made  small  by  suitable  hypotheses,  as  in  [13].  Thus 
we  have  shown 

Theorem  4 . 1 : Assume  that  as  Ax  , Ay  0 , the  discrete 

solution  obtained  from  (4. 3-4. 5)  remains  uniformly  bounded, 

and  converges  in  measure  to  a limit  pair  of  functions  p,q, 

with  p of  bounded  variation  in  x,  uniformly  with  respect 

. °° 

to  y.  Then  for  any  nonnegative  <f>  €'  Cq, 

(4.7)  ( (|p3+|q2  )<|>x-pq<j>y  )dxdy  * 0, 

which  is  the  desired  entropy  inequality. 

Next  we  consider  the  order  of  accuracy  of  the  difference 
scheme.  In  regions  where  the  solution  is  smooth  and  which 
are  away  from  the  sonic  line  (p=0),  each  of  the  difference 
equations  is  clearly  of  second  order  accuracy.  Thus  we 
anticipate  no  trouble  from  the  changing  of  the  form  of  the 
difference  equations  near  the  x-boundaries . 

In  regions  where  the  solution  is  not  smooth,  it  still 
makes  sense  to  determine  to  what  order  of  accuracy  a differ- 
ence scheme  approximates  the  weak  form  of  a given  differ- 
ential equation.  In  the  case  of  linear  hyperbolic  systems, 
this  is  sufficient  to  determine  the  order  of  magnitude  of  the 
error  in  distribution  sense  [13,  14].  To  be  precise,  let 


I 
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q,(x,y)  he  the  function  obtained  by  extension  of  the  discrete 

. ,,  k + 1/2  . . . ^ ^ 

variables  ^+1/2’  ^h  1S  Piecewise  constant,  and  is  our  appr- 

k 

oximation  to  q.  Similarly,  let  p^ , k = 0,1,..., N,  be  a 

function  of  x,  equal  to  the  right  side  of  (4.2),  and  let 

k 

p^  = p^(x,y)  be  the  extension  of  the  p^  with  respect  to 
y.  The  function  p^  is  piecewise  linear  in  x,  piecewise 
constant  in  y,  and  is  our  approximation  to  p. 

Theorem  4.2:  Suppose  that  as  Ax,  Ay  -+0  with  Ay/Ax 
fixed,  the  discrete  solution  , q is  uniformly  bounded 
and  that  p^  is  bounded  variation  in  x,  uniformly  with 
respect  to  y.  Then  for  any  4>  € c”. 


(4.8) 


and 


(4.9) 


(phVqhVdxdyl  * 0(Ax2  + Ay2), 


(qhVPhVdXdyl  * 0 ( Ax2  + Ay2  ) , 


F'roof : Let  us  begin  with  the  second  term  in  (4.8) 


qh^ydxdy  = qh ,y ^dxdy 


=.I  (qk+1/2-qk-1/2) 

j ,k  qj  + l./2  qj  + l/2; 


t ( j +1 ) Ax 


j Ax 


<t>  (x  , (k+ j)Ay  )dx 


. v . k+1/2  k-1/2 . 1 fv 

'j!k  ,i+1/2‘,i+1/2  Jjax 


( j + l)Ax 


(k+1 ) Ay  2 2 

4>(x,y)dxdy  + 0(Ax  +Ay^) 

kAy 
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, k+1/2  k-1/2, 

lx4y  i 


j 

AxAy  l 
j >k 


k+1/2  k-1/2 

!~q- 

'Ay' 

k+1/2  k-1/2 

Lli 

Ay 


(k+1 ) Ay 


kAy 


4>(  (j+j)Ax,y)dy+0(Ax2  + Ay2) 


/ K'+O./  C.  _ -IX  .I-  / v 

kqi+l/2~qj+l/2;  _k 


“1+1/2  + 0(Ax2+Ay2 ) > 


(4.10) 


using  the  piecewise  constant  form  of  and  the  smoothness 

of  <p . We  note  in  passing  that  for  the  linear  relation  (4.9), 
both  terms  are  treated  in  essentially  this  manner  and  (4.7) 
follows  from  the  difference  equation  (4.4).  In  (4.8),  how, 
ever,  the  quadratic  term  is  of  central  importance,  and  it  is 
here  that  the  peculiear  interpretation  (4.2)  of  the  p^ 
is  used.  Let  X be  the  space  of  continuous,  piecewise  lin- 
ear functions  of  x,  with  nodes  at  the  -mesh  points  jAx. 

k 2 

For  each  k,  p^€X.  Now  the  approximation  of  (p  )x  in 

(4.5)  is  just  what  would  be  obtained  (with  respect  to  the 

x-variable)  by  Galerkin's  method  with  the  space  X.  (This 

is  also  an  explanation  of  why  the  entropy  inequality  (4.6) 

• )c  . • • • 

was  obtained. ) Let  y £X  be  the  piecewise  linear  extension 

]< 

of  the  aj  + i/2’  Phen 

. , (k+l)Ay  9 9 

(4.11)  | y (x)  - j—  <j>(x,y)dy|  = 0(Ax  +Ay^) 

Ay  kAy 

uniformly  in  x.  Using  the  bounded  variation  of  p^,  we 


have  from  (4.11) 
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* 


i 


(k+1 ) Ay 
kAy 


(4.12) 


Ph^dxdy  = j|<p£>x  h 

+ CKAx^+Ay^) 


-Ay 


<p  ( x , y )dydx 

k*  2 k . 2.  .2 


■ AxAy  £ 

j 


(yj43/2~Ui-l/2)  (vJj-i-l/2'l'tJi-t-l/2't'Pi-l/2) 
3 Ax 


aj+l/2+0(Ax2+Ay2) 


Finally,  multiply  (4.5)  by  AxAy  aj  + ^/2  anc*  sum  over 

j,k.  Using  the  smoothness  of  4> , the  aritificial  viscosity 

3 3 

term  gives  a continuation  of  0(Ax  +Ay  ),  so  (4.8)  follows 
by  comparison  with  (4.10)  and  (4.12).  This  concludes  the 
proof. 

Finally,  we  describe  a potentially  helpful  modification 
of  the  difference  scheme  (4. 3-4. 5),  corresponding  to  the 
simple  change  of  the  variable  p p+3  for  a constant  3 
to  be  determined  empirically.  The  boundary  conditions  at 
x = 0 and  x = L become  p = 3-Pq,  and  an  additional  term 

k k 

-2ft 

Ax 


is  added  to  the  left  side  of  (4.5).  This  transformation 
is  equivalent,  in  the  difference  scheme,  to  adding  a disper- 
sion  term  of  order  3Ax  Pxxx*  While  such  a term  will  not 
aid  in  the  generation  of  entropy  in  the  vicinity  of  a shock, 
it  may  help  in  controlling  the  variation  of  p^.  This  may  be 
a serious  problem  in  some  applications;  in  particular,  the 
fourth-order  viscosity  invariably  leads  to  oscillatory,  dis- 
crete shock  profiles. 
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Appendix 

On  the  nonexistence  of  weak  shock-type  discontinuities  in 
elliptic  regions. 


Levi  Lustman,  Nima  Geffen 

I 
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Consider  N conservation  equations: 


V0V  = (x,y,u)  + B^(x,y,u)  =0  i = 


U = {Uj } 


j = 1, . . . ,N 


V1  = (A1,  B1!  Cx,y,u) 


Let  A1,  B1  be  continuously  differentiable  (6C1)  in  all 
their  variables  and:  


i . 3A1  Ri  _ 9Ba 
j " “5117  ’ ' aTTT 


eq.  (1)  can  be  written 


a!  |ni  + |h! 

3 9x  ] 9y 


with  summation  convention  on  repeated  indices. 

Let  (Al)  be  elliptic  in  a region  E in  u space 
for  Cx,y)  in  ftj,,  for  which: 


IN 

l ( AA^  - yB^v3  = 0 ■*  v1  = 0 

i = l 3 3 


for  A,  y real  and  not  both  zero. 


I 


I 

I 
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Equation  (A3)  is  a requirement  that  no  real  charac- 
teristics exist,  or  equivalently,  that  (Al)  is  elliptic. 

It  is  symmetric  in  (x,y)  and  N has  to  be  even,  i.e.  N=2m. 
Let : 


ux  € E 

and : 


S(u^)  be  the  points  on  the  shock  polar  emanat- 
ing at  u^  (for  which  a shock  jump  from  u-^  is  compatible 
with  the  conservation  laws  (Al). 


< 

■ 


' 


! 


j 
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Theorem  (Al) 

is  an  isolated  point  of  sCu^)  (hence  no  'weak' 
shocks  near  u-^€E  are  possible). 

Note  that  by  continuity,  there  exists  a neighborhood  K 
of  u1  such  that  for  all  ( Z ,m) , (u  ,v)  ( l jiiOeK^A1  ,B1(u)  ( £ ,m) 
will  leave  (Al)  elliptic.  Let  K be  con- 
vex, the  claim  is  that  there  is  no  jump  connection  between 
u^  and  uk  € K for  k i 1. 

Proof : 

Assume:  u2  e SCu-^),  u2  € K,  0 i 0. 

Then : 


B1(x,y,u_1)  - B1(x,y,u2) 
A1 (x  ,y ,u^)  - A1(x,y,u2) 


Let,  for  0*0*1: 

A1(6)  = A1(x,y,(l-0)u1+  u 2 ) 

B^(0)  = B^(x,y , (1-0)11^  + 0u2) 
we  get : 

(A1,B1)(0)  = (A1,Bi)(x,y,u1) 
(A1,B1)(1)  = (A1,B1)(x,y,u2) 


I 
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and,  using  the  mean-value  theorem,  there  exist 
®A’  ®B  = l»***»lj  such  that: 


Al^,y,u_2)  - A^'Cx^jU^)  = 


dA 


0 = 0; 


B Cx ,y ,u2 ) - B Cx,y,u^)  = ^ — 


0 = 0; 


But  we  have : 

, i N 


dA 

d0 


= I A^  (x,y,Cl-0)u,  + 0u  ) (u9  - 
i=l  3 ~2  “2 


-l5 


i N 


dB  i 

- l B-i  (x,y , ci-0)u,  + 0U9)  (u9  - 
j = l J z 


Hi) 


and  since: 


(1-0 )ux  + 0u2  £ K 


we  have : 


for: 


0 ■ j,  (4A) 

= Aj  (x,y,(l-0A)u1  + 0Au2> 

= Bj  (x,y,(l-0B)u1  + 6bu2) 


i = 1,. . . ,N 
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letting : 


we  get : 


v = u 2 - u i £ K 

€ K 


contradicting  the  ellipticity  condition  (A3). 

This  proves  the  non  existence  of  arbitrarily  weak  shocks 
in  elliptic  fields,  precludes  the  possibility  of  closed  shock 
polars  in  u space  (e.g.  the  hodograph  plane)  and  of  "shock- 
ed" elliptic  flows  about  obstacles  (with  shock  strength  van- 
ishing at  infinite). 


